Analysis of Sudan survey data 2012 - 2013

This R-markdown document contains all r-code used to carry out the analysis for the paper on fish catch rates and species densities (species richness) along the Red Sea coast of Sudan. The analysis is based on three surveys: November 2012, May 2013 and November 2013.

Here the analysis is structured around the 7 management regions of the Sudanese Red Sea coast as follows:

All code and data are stored on GitHub

Libraries, data etc.

Libraries, dependencies and functions

Load and manipulate data

Reading catch, station and traits data.

Neither station data, nor catch data has a complete depth record, but by combining depth para from each we can get a complete depth parameter for all stations.

Sudan map and management areas

Loads the map data for Sudan, loads the management areas from shape-file etc.

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/eriko/Documents/GitHub/Sudan2019/Sudan-master/sudan_management_areas", layer: "sudan_regions"
## with 7 features
## It has 3 fields
## Integer64 fields read as strings:  id

Allocating catch positions to management areas

Adds ‘Area’ to each line in the catch table.

Modifying dataset and creating subset for stations with catches

Adding traits to ‘catch’ data set (without 0-catch stations)

Some more data wrangling that adds the traits from the traits table to the catch data, using only the station with catches (there are no traits for ‘NOCATCH’ species).

Also, only select traps, gillnets and handlines as these were the only gear with sufficient numbers and consistent use to be analyzed.

Lastly adds number of gear deployed at each station to each line.

Bathymetric map of Sudan

Fig. 1 in MS

A nice, bathymetric map of Sudan with management areas overlaid.

Currently fonts and colours are adjusted for use on the ICES ASC 2019 poster, removing the bathy-contour lines (se ‘lwd’ to 0), and varying the colours of management area text to white for areas 1-4.

## quartz_off_screen 
##                 2

Stations plotted on maps, pr survey

Faceted map plotting position of all catch stations for each of the three surveys.

MApoints<-data.frame(lon=c(37.2, 37.53, 37.5, 37.4, 37.3, 37.9, 38.2), lat=c(21.7, 20.92, 20.3, 19.9, 19, 19.3, 18.5), area=c(1,2,3,4,5,6,7))

sudan.map<-ggplot(catch, aes(x=lon, y=lat, group=survey_m)) + geom_point(shape="+", size=6, aes(group=survey_m, colour= survey_m)) + facet_wrap(~survey_m, ncol=3) + geom_path(data=ManageAreas.f,colour="dodgerblue3", aes(x=long, y=lat, group=id)) + geom_polygon(data=world, colour="gray35", fill="gray85", aes(x=long, y=lat, group=survey_m)) +  coord_cartesian(xlim = c(36.5, 39), ylim=c(17.7, 22.5)) + theme_classic() + geom_point(data=ps, size=2, colour="gray35", aes(x=x, y=y)) + geom_text(data=ps, label="PZU", size=3,  hjust=1, vjust=-1.2,  aes(x=x, y=y, group=survey_m)) +geom_text(data=MApoints, aes(x=lon, y=lat, group=area, label=area)) +scale_colour_brewer(type = "seq", palette = "Dark2", name="Survey month & year", labels=c("Nov. 2012", "May. 2013", "Nov. 2013")) + theme(legend.position="none")
sudan.map

ggsave("./figs/station_map.tiff")
## Saving 7 x 5 in image

Station information

Station table

Table describing the sampling effort in each area pr survey, number of different gear types, max, min and average depth, number of traps with / without catches.

survey id Ntraps Nhl NGn TBhrs HLhrs GNhrs DepthAvg DepthMax DepthMin
2012901 1 22 0 0 694.066 0.000 0.000 42.45455 142 13
2012901 2 54 3 3 721.549 78.000 62.000 40.89744 71 0
2012901 3 26 0 1 451.183 0.000 14.000 30.80769 95 8
2012901 4 5 0 0 77.084 0.000 0.000 22.60000 54 10
2012901 5 31 0 4 677.896 0.000 78.317 21.32258 30 7
2012901 6 36 0 1 850.417 0.000 38.250 32.38889 88 0
2012901 7 31 0 8 712.200 0.000 135.849 31.06667 66 5
2013002 1 29 0 1 420.163 0.000 24.000 31.48148 70 5
2013002 2 81 3 5 1102.379 377.500 221.569 27.09091 145 5
2013002 3 32 0 1 545.811 0.000 24.000 28.87500 60 0
2013002 4 13 0 10 160.266 0.000 137.815 29.50000 67 9
2013002 5 33 2 5 208.899 192.000 195.766 20.46154 50 7
2013002 6 45 1 2 642.413 156.000 78.000 34.70270 88 9
2013002 7 39 2 0 666.410 186.000 0.000 33.46154 76 5
2013005 1 23 1 4 271.551 3.000 84.000 29.82353 80 10
2013005 2 57 2 6 500.101 111.000 156.000 38.27273 80 7
2013005 3 9 0 2 123.099 0.000 36.000 26.50000 70 9
2013005 4 16 2 4 151.184 4.500 142.767 32.90000 68 12
2013005 5 30 2 3 317.498 7.000 146.947 25.52381 65 11
2013005 6 40 4 2 443.983 31.501 71.915 40.14815 89 6
2013005 7 22 0 2 171.784 0.000 203.171 34.75000 54 11

Species table

Number of fish (organized by family and species) caught by Gillnet or Traps for each survey, and in total across all surveys and gears.

  • (new table for revision 2020) *
    fam_name Sci_name 2012901_GN 2012901_TB 2013002_GN 2013002_TB 2013005_GN 2013005_TB sum
    ACANTHURIDAE Acanthurus nigricans 0 10 0 24 6 49 89
    ACANTHURIDAE Acanthurus nigricauda 0 0 0 12 0 0 12
    SOLEIDAE Achirus sp. * 0 0 0 0 1 0 1
    HOLOCENTRIDAE Adioryx ruber 0 0 0 11 0 4 15
    HOLOCENTRIDAE Adioryx spinifer 0 0 0 0 1 0 1
    SERRANIDAE Aethaloperca rogaa 0 2 0 10 0 1 13
    ALBULIDAE Albula vulpes 0 0 0 0 5 0 5
    CARANGIDAE Alectis indicus 0 0 0 0 1 0 1
    CARANGIDAE Alepes vari 0 0 0 0 4 0 4
    SPARIDAE Argyrops filamentosus 0 0 0 7 0 0 7
    SPARIDAE Argyrops sp. 0 20 0 0 0 0 20
    SPARIDAE Argyrops spinifer 0 0 0 6 0 5 11
    ARIIDAE Arius heudelotii 0 6 0 0 0 0 6
    ARIIDAE Arius thalassinus 0 0 0 2 1 0 3
    SCOMBRIDAE Auxis thazard 14 0 0 0 0 0 14
    BALISTIDAE Balistapus undulatus 0 0 0 0 0 1 1
    BALISTIDAE Balistoides viridescens 0 0 0 1 0 0 1
    BOTHIDAE Bothus pantherinus 0 0 0 0 1 0 1
    CAESIONIDAE Caesio caerulaurea 0 0 15 0 0 0 15
    CAESIONIDAE Caesio suevicus 0 0 0 0 1 0 1
    CARANGIDAE Carangoides armatus 0 0 0 0 4 0 4
    CARANGIDAE Carangoides bajad 13 1 6 14 49 0 83
    CARANGIDAE Carangoides ferdau 0 0 0 1 7 0 8
    CARANGIDAE Carangoides fulvoguttatus 0 0 0 0 16 0 16
    CARANGIDAE Carangoides sp. 0 0 0 0 1 0 1
    CARANGIDAE Caranx (Caranx) melampygus 2 0 0 1 6 0 9
    CARANGIDAE Caranx (Caranx) sexfasciatus 16 0 22 2 43 0 83
    CARANGIDAE Caranx (Gnathanodon) speciosus 0 0 0 0 4 0 4
    CARANGIDAE Caranx ignobilis 0 0 0 1 1 0 2
    CARANGIDAE Caranx sp. 0 0 0 0 1 0 1
    CARCHARHINIDAE Carcharhinus albimarginatus 2 0 0 0 0 0 2
    S H A R K S Carcharhinus melanopterus 9 0 0 0 5 0 14
    S H A R K S Carcharhinus wheeleri 0 0 1 0 1 0 2
    SERRANIDAE Cephalopholis argus 0 1 0 0 0 0 1
    SERRANIDAE Cephalopholis miniatus ? 0 0 0 0 1 0 1
    CHAETODONTIDAE Chaetodon auriga 0 3 0 0 0 0 3
    CHAETODONTIDAE Chaetodon semilarvatus 0 0 16 2 0 0 18
    CHANIDAE Chanos chanos 0 0 0 0 1 0 1
    LABRIDAE Cheilinus lunulatus 0 0 0 1 0 0 1
    LABRIDAE Cheilinus quinquecintus 0 1 0 0 0 0 1
    CHIROCENTRIDAE Chirocentrus dorab 15 0 17 0 60 0 92
    CARANGIDAE Chorinemus lysan 0 0 0 6 0 0 6
    PLATYCEPHALIDAE Cociella crocodila 0 0 0 0 1 0 1
    MUGILIDAE Crenimugil crenilabis 0 0 1 0 0 0 1
    CARANGIDAE Decapterus macarellus 0 0 1 0 0 0 1
    SCOMBRIDAE Decapterus russelli 0 0 1 4 0 0 5
    HAEMULIDAE Diagramma pictum 0 1 0 0 0 0 1
    DIODONTIDAE Diodon hystrix 0 0 0 0 3 0 3
    ECHENEIDIDAE Echeneis naucrates 0 0 0 3 2 0 5
    SCOMBRIDAE Elagatis bipinnulata 0 0 6 0 0 0 6
    SERRANIDAE Epinephelus chlorostigma 0 0 0 1 0 0 1
    SERRANIDAE Epinephelus fuscoguttatus 0 7 0 7 1 2 17
    SERRANIDAE Epinephelus sexfasciatus 0 0 0 1 0 1 2
    SERRANIDAE Epinephelus summana 0 0 2 1 0 0 3
    SERRANIDAE Epinephelus tauvina 1 11 4 4 3 2 25
    SCOMBRIDAE Euthynnus affinis 0 0 4 0 3 0 7
    JUVENILES FISHJUV 9 0 0 1 0 0 10
    GERREIDAE Gerres oyena 0 0 1 0 6 0 7
    SCOMBRIDAE Grammatorcynus bicarinatus 0 0 3 0 0 0 3
    SCOMBRIDAE Grammatorcynus bilineatus 26 2 2 0 5 0 35
    LETHRINIDAE Gymnocranius microdon 0 0 1 0 0 0 1
    LETHRINIDAE Gymnocranius robinsoni 0 0 0 2 3 0 5
    SCOMBRIDAE Gymnosarda unicolor 0 0 0 0 10 0 10
    MURAENIDAE Gymnothorax flavimarginatus 0 3 0 0 0 0 3
    MURAENIDAE Gymnothorax javanicus 0 28 0 11 0 1 40
    HAEMULIDAE Hectromynicus pictus 0 0 0 4 0 0 4
    HEMIRAMPHIDAE Hemiramphus far 0 0 0 0 1 0 1
    SCARIDAE Hipposcarus harid 0 0 4 0 4 0 8
    SCOMBRIDAE Katsuwonus pelamis 1 0 0 0 0 0 1
    KYPHOSIDAE Kyphosus vaigiensis 0 0 0 0 10 0 10
    LETHRINIDAE Lethrinus elongatus * 0 10 1 17 7 5 40
    LETHRINIDAE Lethrinus harak 0 0 0 0 10 0 10
    LETHRINIDAE Lethrinus lentjan 3 33 6 47 31 28 148
    LETHRINIDAE Lethrinus mahsena 0 16 2 19 0 18 55
    LETHRINIDAE Lethrinus mahsenoides * 0 4 0 0 6 0 10
    LETHRINIDAE Lethrinus nebulosus 0 0 0 0 0 1 1
    LETHRINIDAE Lethrinus ramak 0 0 1 1 0 1 3
    LETHRINIDAE Lethrinus xanthochilus 0 0 0 2 0 1 3
    LUTJANIDAE Lutjanus argentimaculatus 0 0 0 1 0 1 2
    LUTJANIDAE Lutjanus bohar 0 32 10 69 3 8 122
    LUTJANIDAE Lutjanus cf fluviflamma 0 0 0 0 2 0 2
    LUTJANIDAE Lutjanus coccineus * 0 0 0 0 0 1 1
    LUTJANIDAE Lutjanus ehrenbergii 6 0 11 0 7 0 24
    LUTJANIDAE Lutjanus fulviflamma 0 0 0 0 2 0 2
    LUTJANIDAE Lutjanus gibbus 0 38 0 51 1 28 118
    LUTJANIDAE Lutjanus kasmira 0 4 0 3 0 2 9
    LUTJANIDAE Lutjanus monostigma 2 4 0 3 0 1 10
    LUTJANIDAE Lutjanus rivulatus 0 0 0 1 0 0 1
    LUTJANIDAE Lutjanus sebae 0 0 0 1 0 0 1
    LUTJANIDAE Lutjanus sp. 0 1 0 0 0 0 1
    LUTJANIDAE Macolor niger 0 0 0 2 1 0 3
    LETHRINIDAE Monotaxis grandoculis 0 0 1 0 0 0 1
    MULLIDAE Mulloides flavolineatus 0 0 1 0 0 0 1
    MULLIDAE Mulloides vanicolensis 0 0 0 0 1 0 1
    HOLOCENTRIDAE Myripristis murdjan 0 0 3 0 4 0 7
    ACANTHURIDAE Naso hexacanthus 0 0 18 18 37 0 73
    ACANTHURIDAE Naso lituratus 0 0 0 0 1 0 1
    NO CATCH NO CATCH 3 0 0 0 0 0 3
    LUTJANIDAE Paracaesio sordius 0 2 0 0 0 0 2
    EPHIPPIDAE Platax boersi 0 1 0 0 0 0 1
    EPHIPPIDAE Platax orbicularis 0 2 0 2 1 2 7
    HAEMULIDAE Plectorhinchus gaterinus 0 7 1 0 3 0 11
    POMADASYIDAE (HAEMULIDAE) Plectorhinchus pictus 0 0 0 0 1 0 1
    HAEMULIDAE Plectorhynchus pictus 0 0 0 0 1 0 1
    HAEMULIDAE Plectorhynchus schotaf 0 0 0 0 1 0 1
    SERRANIDAE Plectropomus pessuliferus 0 2 0 2 0 0 4
    PRIACANTHIDAE Priacanthus hamrur 0 0 4 0 0 0 4
    LUTJANIDAE Pristipomoides multidens 0 9 0 7 0 0 16
    BALISTIDAE Pseudobalistes flavimarginatus 0 1 0 0 0 0 1
    SCOMBRIDAE Rastrelliger kanagurta 0 0 6 0 21 0 27
    SCOMBRIDAE Sarda orientalis 0 0 1 0 0 0 1
    HOLOCENTRIDAE Sargocentron spiniferum 0 13 1 13 0 5 32
    SCARIDAE Scarus ferrugineus 0 0 1 0 0 0 1
    SCARIDAE Scarus frenatus 0 0 2 0 2 0 4
    SCARIDAE Scarus ghobban 0 0 0 0 1 0 1
    SCOMBRIDAE Scomber japonicus 0 0 0 1 0 0 1
    CARANGIDAE Scomberoides lysan 15 20 60 18 45 0 158
    CARANGIDAE Scomberoides tol 0 0 4 0 74 0 78
    SCOMBRIDAE Scomberomorus commerson 35 0 0 0 0 0 35
    SCOMBRIDAE Scomberomorus lineolatus 0 0 0 0 9 0 9
    SCOMBRIDAE Scomberomorus tritor 2 0 0 0 0 0 2
    SIGANIDAE Siganus argenteus 0 0 3 0 0 0 3
    SIGANIDAE Siganus luridus 1 0 1 0 0 0 2
    SIGANIDAE Siganus rivulatus 0 0 0 0 1 0 1
    SIGANIDAE Siganus stellatus 0 0 1 0 0 2 3
    SPARIDAE Sparus sp. 0 0 0 6 0 0 6
    SPHYRAENIDAE Sphyraena forsteri 0 0 1 0 0 0 1
    SPHYRAENIDAE Sphyraena jello 3 0 0 0 0 0 3
    SPHYRAENIDAE Sphyraena putnamie 0 0 0 0 1 0 1
    SPHYRAENIDAE Sphyraena qenie 0 0 6 0 4 0 10
    CARCHARHINIDAE Sphyrna lewini 1 0 0 0 0 0 1
    R A Y S Taeniura lymma 0 0 0 0 1 0 1
    SCOMBRIDAE Thunnus albacares 11 0 0 0 0 0 11
    CARCHARHINIDAE Triaenodon obesus 0 3 0 10 0 0 13
    BELONIDAE Tylosurus choram 2 0 4 0 0 0 6
    MUGILIDAE Valamugil engeli 0 0 0 0 1 0 1
    SERRANIDAE Variola louti 0 0 0 2 0 0 2

Number of set traps per depth

Violin plot of depths at wich traps were set. I prefer this visualization of depth ranges of the traps, rather than a (boring) table.

c2 <- catch %>% dplyr::group_by(ss) %>% dplyr::summarise(first(depth), first(gear), first(id), first(survey), first(survey_m))
colnames(c2) <- c("ss", "depth", "gear", "id", "survey", "survey_m")

dp <- ggplot(subset(c2, gear=="TB" & depth!=0), aes(as.factor(id), depth))

dp <- dp +  geom_violin(scale="width", aes(fill=as.factor(id)), draw_quantiles = c(0.5)) + scale_y_reverse() + labs(x="Management areas", y="Depth (m)") + scale_fill_manual(name="Management area",values=sequential_hcl(n = 7, h = c(0, -100), c = c(80, NA, 40), l = c(53, 75), power = c(1, 1), rev = TRUE, register = ), guide=FALSE ) +theme_bw()  +facet_wrap(~as.factor(survey_m))

dp

tiff(file="./figs/trap_area_depth.tiff", width=3000, height=1500, res=400, pointsize=10, compression=c("none"))
dp
dev.off()
## quartz_off_screen 
##                 2

Test differneces in depths in each area between years and between areas

## 
##  Shapiro-Wilk normality test
## 
## data:  c2$depth
## W = 0.85976, p-value < 2.2e-16

## 
##  Kruskal-Wallis rank sum test
## 
## data:  c2$depth and c2$id
## Kruskal-Wallis chi-squared = 33.554, df = 6, p-value = 8.201e-06
## Warning in posthoc.kruskal.nemenyi.test.default(c2$depth, c2$id, "Chisq"): Ties
## are present. Chi-sq was corrected for ties.
## 
##  Pairwise comparisons using Nemenyi-test with Chi-squared    
##                        approximation for independent samples 
## 
## data:  c2$depth and c2$id 
## 
##   1       2       3       4       5       6      
## 2 0.99215 -       -       -       -       -      
## 3 0.59634 0.82901 -       -       -       -      
## 4 0.70243 0.89758 1.00000 -       -       -      
## 5 0.00097 0.00055 0.51581 0.71502 -       -      
## 6 0.93649 0.99855 0.98040 0.98709 0.02090 -      
## 7 0.99461 1.00000 0.90318 0.93635 0.00700 0.99962
## 
## P value adjustment method: none
## [1] "2012901"
## 
##  Kruskal-Wallis rank sum test
## 
## data:  ci$depth and ci$id
## Kruskal-Wallis chi-squared = 24.982, df = 6, p-value = 0.0003441
## Warning in posthoc.kruskal.nemenyi.test.default(ci$depth, ci$id, "Chisq"): Ties
## are present. Chi-sq was corrected for ties.
## 
##  Pairwise comparisons using Nemenyi-test with Chi-squared    
##                        approximation for independent samples 
## 
## data:  ci$depth and ci$id 
## 
##   1     2     3     4     5     6    
## 2 1.000 -     -     -     -     -    
## 3 0.510 0.560 -     -     -     -    
## 4 0.604 0.689 0.992 -     -     -    
## 5 0.035 0.016 0.953 1.000 -     -    
## 6 0.386 0.384 1.000 0.993 0.941 -    
## 7 0.842 0.909 0.991 0.931 0.429 0.979
## 
## P value adjustment method: none 
## [1] "2013002"
## 
##  Kruskal-Wallis rank sum test
## 
## data:  ci$depth and ci$id
## Kruskal-Wallis chi-squared = 14.173, df = 6, p-value = 0.02776
## Warning in posthoc.kruskal.nemenyi.test.default(ci$depth, ci$id, "Chisq"): Ties
## are present. Chi-sq was corrected for ties.
## 
##  Pairwise comparisons using Nemenyi-test with Chi-squared    
##                        approximation for independent samples 
## 
## data:  ci$depth and ci$id 
## 
##   1     2     3     4     5     6    
## 2 0.983 -     -     -     -     -    
## 3 0.998 1.000 -     -     -     -    
## 4 1.000 1.000 1.000 -     -     -    
## 5 0.377 0.637 0.741 0.740 -     -    
## 6 1.000 0.781 0.960 0.992 0.088 -    
## 7 1.000 0.917 0.989 0.998 0.186 1.000
## 
## P value adjustment method: none 
## [1] "2013005"
## 
##  Kruskal-Wallis rank sum test
## 
## data:  ci$depth and ci$id
## Kruskal-Wallis chi-squared = 9.081, df = 6, p-value = 0.1691
## Warning in posthoc.kruskal.nemenyi.test.default(ci$depth, ci$id, "Chisq"): Ties
## are present. Chi-sq was corrected for ties.
## 
##  Pairwise comparisons using Nemenyi-test with Chi-squared    
##                        approximation for independent samples 
## 
## data:  ci$depth and ci$id 
## 
##   1    2    3    4    5    6   
## 2 1.00 -    -    -    -    -   
## 3 0.99 0.95 -    -    -    -   
## 4 0.99 0.94 1.00 -    -    -   
## 5 0.66 0.28 1.00 0.99 -    -   
## 6 1.00 0.99 1.00 1.00 0.77 -   
## 7 1.00 1.00 0.99 1.00 0.79 1.00
## 
## P value adjustment method: none

Depth data were normally distributed and hence the difference in depth could be analysed with ANOVA.

Of the 63 survey - area comparisons, only two depth distributions were significantly different (Kruskal-Wallis rank-sum test and Nemenyi post-hoc test, p<0.05); in the Nov. 2012 between areas 5:1 and 5:2, - traps in areas 1 and 2 were set significantly deeper than traps in area 5. However, the K-W test also showed significant differences for the May 2013 survey, but the Nemenyi test could not detect any significant differences between the areas.

Duration of trap sets

Mean, min & max hrs. of fishing

id maxhrs minhrs meanhrs
1 46.417 11.500 20.48641
2 44.500 7.917 16.47669
3 24.250 8.167 16.86333
4 27.033 13.033 16.11076
5 40.200 11.667 17.79179
6 30.083 8.000 19.00436
7 31.866 10.417 19.47891 
##        ss                  id             fhrs       
##  Min.   :2.013e+09   Min.   :1.000   Min.   : 7.917  
##  1st Qu.:2.013e+09   1st Qu.:2.000   1st Qu.:15.425  
##  Median :2.013e+09   Median :4.000   Median :16.150  
##  Mean   :2.013e+09   Mean   :3.909   Mean   :17.984  
##  3rd Qu.:2.013e+09   3rd Qu.:6.000   3rd Qu.:18.700  
##  Max.   :2.013e+09   Max.   :7.000   Max.   :46.417

Fishing hours (traps) by survey and area

Boxplot of median with mean values plotted as red circle. The average trap soak-time (fishing hours) across all surveys and areas was 17.984 hrs with a median of 16.15, but the November 2012 survey had more varied and higher soak times than the May 2013 and Nov. 2013.

Test differences in soak times (traps)

## 
##  Shapiro-Wilk normality test
## 
## data:  c3$Fhrs
## W = 0.67297, p-value < 2.2e-16

## 
##  Kruskal-Wallis rank sum test
## 
## data:  c3$Fhrs and c3$id
## Kruskal-Wallis chi-squared = 68.971, df = 6, p-value = 6.645e-13
## Warning in posthoc.kruskal.nemenyi.test.default(c3$Fhrs, c3$id, "Chisq"): Ties
## are present. Chi-sq was corrected for ties.
## 
##  Pairwise comparisons using Nemenyi-test with Chi-squared    
##                        approximation for independent samples 
## 
## data:  c3$Fhrs and c3$id 
## 
##   1       2       3       4       5       6      
## 2 0.89835 -       -       -       -       -      
## 3 1.00000 0.88439 -       -       -       -      
## 4 0.61831 0.95548 0.60026 -       -       -      
## 5 0.99997 0.65429 1.00000 0.42582 -       -      
## 6 0.14066 1.3e-05 0.20277 0.00167 0.19214 -      
## 7 0.01625 2.0e-07 0.02891 0.00012 0.02234 0.97141
## 
## P value adjustment method: none
## [1] "2012901"
## 
##  Kruskal-Wallis rank sum test
## 
## data:  ci$Fhrs and ci$id
## Kruskal-Wallis chi-squared = 60.806, df = 6, p-value = 3.087e-11
## Warning in posthoc.kruskal.nemenyi.test.default(ci$Fhrs, ci$id, "Chisq"): Ties
## are present. Chi-sq was corrected for ties.
## 
##  Pairwise comparisons using Nemenyi-test with Chi-squared    
##                        approximation for independent samples 
## 
## data:  ci$Fhrs and ci$id 
## 
##   1       2       3       4       5       6      
## 2 0.07376 -       -       -       -       -      
## 3 0.06785 0.99940 -       -       -       -      
## 4 0.18427 0.94372 0.98714 -       -       -      
## 5 0.96585 0.48177 0.39946 0.46876 -       -      
## 6 0.99285 0.00023 0.00085 0.04610 0.49836 -      
## 7 0.96104 8.6e-05 0.00033 0.02854 0.32747 0.99990
## 
## P value adjustment method: none 
## [1] "2013002"
## 
##  Kruskal-Wallis rank sum test
## 
## data:  ci$Fhrs and ci$id
## Kruskal-Wallis chi-squared = 42.845, df = 6, p-value = 1.252e-07
## Warning in posthoc.kruskal.nemenyi.test.default(ci$Fhrs, ci$id, "Chisq"): Ties
## are present. Chi-sq was corrected for ties.
## 
##  Pairwise comparisons using Nemenyi-test with Chi-squared    
##                        approximation for independent samples 
## 
## data:  ci$Fhrs and ci$id 
## 
##   1       2       3       4       5       6      
## 2 0.99993 -       -       -       -       -      
## 3 0.23422 0.14820 -       -       -       -      
## 4 0.99997 1.00000 0.74873 -       -       -      
## 5 0.37056 0.28089 0.99998 0.85136 -       -      
## 6 0.00475 0.00024 0.94518 0.20902 0.83861 -      
## 7 0.06190 0.01735 0.99966 0.51167 0.99399 0.99505
## 
## P value adjustment method: none 
## [1] "2013005"
## 
##  Kruskal-Wallis rank sum test
## 
## data:  ci$Fhrs and ci$id
## Kruskal-Wallis chi-squared = 27.538, df = 6, p-value = 0.0001148
## Warning in posthoc.kruskal.nemenyi.test.default(ci$Fhrs, ci$id, "Chisq"): Ties
## are present. Chi-sq was corrected for ties.
## 
##  Pairwise comparisons using Nemenyi-test with Chi-squared    
##                        approximation for independent samples 
## 
## data:  ci$Fhrs and ci$id 
## 
##   1     2     3     4     5     6    
## 2 0.375 -     -     -     -     -    
## 3 0.307 0.966 -     -     -     -    
## 4 0.962 0.995 0.870 -     -     -    
## 5 0.302 1.000 0.995 0.973 -     -    
## 6 0.979 0.830 0.627 1.000 0.719 -    
## 7 0.958 0.015 0.045 0.496 0.015 0.426
## 
## P value adjustment method: none

Of the 63 area- survey combinations 12 were significantly different (Kruskal Wallis rank sum test with Nemenyi post-hoc test, p<0.05). For the Nov 2012 survey soak time in areas 6 and 7 was significantly longer than in areas 2, 3 and 4. For the May 2013 survey soak times in area 6 were significantly longer than in areas 1 and 2, while soak time in area 7 was significantly longer than in area 2. For the Nov 2013 survey soak time in area 7 was significantly longer than in areas 2, 3 and 5.

Analysis of catches by traits and other factors

Catch rates pr area & survey

ca.plot<-ggplot(subset(catch, gear=="TB" & depth!=0), aes(as.factor(id), CPUEw))

ca.plot <- ca.plot + geom_violin(scale="width", aes(fill=as.factor(id)), draw_quantiles = c(0.5))  + labs(x="Management areas", y="CPUE kg/hrs fishing") + scale_fill_manual(name="Management area",values=sequential_hcl(n = 7, h = c(0, -100), c = c(80, NA, 40), l = c(53, 75), power = c(1, 1), rev = TRUE, register = ), guide=FALSE ) +theme_bw() +facet_wrap(~as.factor(survey_m))
ca.plot

tiff(file="./figs/catch_rates_traps_area_survey.tiff", width=3000, height=2000, res=400, pointsize=10, compression=c("none"))
ca.plot
dev.off()
## quartz_off_screen 
##                 2

CPUE by gear type - by family and trophic group

# by genus (families)
catch.3gear<-subset(catch.unmod, gear=="GN" | gear=="HL" | gear =="TB")

gear.plot <- ggplot(catch.3gear, aes(capwords(tolower(as.character(catch.3gear$FamGroup))), fill=survey_m)) + geom_bar(binwidt=0.2, aes( weight=weight)) +theme_bw() + facet_grid(vars(gear), vars(survey_m)) + coord_flip() +xlab("") + ylab("kg")  + guides(fill=guide_legend(title="Survey")) + theme(legend.position = "none")
## Warning: Ignoring unknown parameters: binwidt
gear.plot

ggsave("./figs/gearplot_species.tiff")
## Saving 7 x 5 in image
# by trophic group
gear.plot2 <- ggplot(subset(catch.traits, TGShort!="NA"), aes(TGShort, fill=survey_m)) + geom_bar(binwidt=0.2, aes( weight=weight)) +theme_bw() + facet_grid(vars(gear), vars(survey_m)) + coord_flip() +xlab("") + ylab("kg")  + guides(fill=guide_legend(title="Survey")) + theme(legend.position = "none")
## Warning: Ignoring unknown parameters: binwidt
gear.plot2

ggsave("./figs/gearplot_trophic_group.tiff")
## Saving 7 x 5 in image

CPUE by trophic group by area & by survey (year)

FIG 5 in Feb 2019 version of MS

Presents the CPUE pr hr, for each management area and survey, standardized by the number of traps set in each area X survey combination.

Trap catches were dominated by carnivores, occuring in all area X survey combinations, followed by invertivores who also occurred everywhere, albeight in lower densities than the carnivores.

Area 2, 5 and 6 seem to have the highest diversity of trophic groups, but this has to be calculated as functional diversity, using the FD package (or similar). 

trophic.plot <- ggplot(subset(catch.traits, TGShort!="NA" & gear=="TB"), aes(TGShort, CPUEw/NArea))

trophic.plot <- trophic.plot + geom_col(aes(fill=TGShort)) + facet_grid(vars(id), vars(survey_m)) + scale_fill_manual(name="Trophic group",values=qualitative_hcl(n = 5, h = c(0, 295), c = 80, l = 60, register = ) ) + labs(x=" ", y="CPUE (kg/hrs fishing) / (number of traps pr area)") + theme_bw() 
trophic.plot

tiff(file="./figs/trophic_group_cpue.tiff", width=3000, height=3000, res=400, pointsize=10, compression=c("none"))
trophic.plot
dev.off()
## quartz_off_screen 
##                 2

Trophic group CPUE by depth

A modified version of Figure 6 in Feb 2019 MS.

Basically a more complex version of the previous plot. Adds the depth dimension to the plot, and the surveys as colours.

Shows uniform catch rates of carnivores and invertivores from 0-50m.

All surveys seem to overlap, so no difference in trophich group depth dependent catch rates between the surveys.

depth.plot <- ggplot(subset(catch.traits, gear=="TB" & TGShort!="NA"), aes(depth, CPUEw)) 

depth.plot <- depth.plot + geom_point(aes(colour=survey_m))  + facet_grid(vars(id),vars(TGShort)) +  scale_colour_manual(name="Survey",values=qualitative_hcl(n = 3, h = c(0, 295), c = 80, l = 60, register = ) ) + labs(x="Depth (m)", y="CPUE (kg(hrs)")  +theme_bw()  

depth.plot

tiff(file="./figs/trophic_group_cpue_depth.tiff", width=3000, height=3000, res=400, pointsize=10, compression=c("none"))
depth.plot
dev.off()
## quartz_off_screen 
##                 2

Trophic level in Trap catches

TMean <- subset(catch.traits, gear=="TB" & TrophicLevel!="NA") %>% group_by(id, survey_m) %>% summarise(TM=mean(na.omit(TrophicLevel)))
CMean <- subset(catch.traits, gear=="TB" & TrophicLevel!="NA") %>% group_by(id, survey_m) %>% summarise(CM=mean(CPUEw))
TMean$CM <- CMean$CM

t0.plot <- ggplot(subset(catch.traits, gear=="TB" & TrophicLevel!="NA"), aes(TrophicLevel, CPUEw)) 

t0.plot <- t0.plot + geom_point(aes(colour=TGShort)) + geom_vline(data = TMean, aes(xintercept = TM)) + facet_grid(vars(id),vars(survey_m)) +  scale_colour_manual(name="Trophic group",values=qualitative_hcl(n = 5, h = c(0, 295), c = 80, l = 60, register = ) ) + labs(x="Trophic level", y="CPUE (kg(hrs)")  +theme_bw()  
t0.plot

tiff(file="./figs/trophic_level.tiff", width=3000, height=3000, res=400, pointsize=10, compression=c("none"))
t0.plot
dev.off()
## quartz_off_screen 
##                 2

Maxlength in catches

t0.plot <- ggplot(subset(catch.traits, gear=="TB" & TrophicLevel!="NA"), aes(MaxLength, CPUEw)) 

t0.plot <- t0.plot + geom_point(aes(colour=TGShort)) + facet_grid(vars(id),vars(survey_m)) +  scale_colour_manual(name="Trophic group",values=qualitative_hcl(n = 5, h = c(0, 295), c = 80, l = 60, register = ) ) + labs(x="Max lenght(cm)", y="CPUE (kg(hrs)")  +theme_bw()  
t0.plot

tiff(file="./figs/max_length.tiff", width=3000, height=3000, res=400, pointsize=10, compression=c("none"))
t0.plot
dev.off()
## quartz_off_screen 
##                 2

Place in water column

catch.traits$WaterCol[catch.traits$WaterCol=="Demersal"] <- "demersal"
tl.plot <- ggplot(subset(catch.traits, WaterCol!="NA" & gear=="TB"), aes(WaterCol, CPUEw/NArea))

tl.plot <- tl.plot + geom_col(aes(fill=WaterCol)) + facet_grid(vars(id), vars(survey_m)) + scale_fill_manual(name="Place in water column",values=qualitative_hcl(n = 6, h = c(0, 295), c = 80, l = 60, register = ) ) + labs(x=" ", y="CPUE (kg/hrs fishing) / number of traps pr area") + theme_bw() + 
  theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank())
tl.plot

tiff(file="./figs/water_col.tiff", width=3000, height=3000, res=400, pointsize=10, compression=c("none"))
tl.plot
dev.off()
## quartz_off_screen 
##                 2

Diel activity

t2.plot <- ggplot(subset(catch.traits, DielActivity!="#N/A" & gear=="TB"), aes(DielActivity, CPUEw/NArea))

t2.plot <- t2.plot + geom_col(aes(fill=DielActivity)) + facet_grid(vars(id), vars(survey_m)) + scale_fill_manual(name="Diel activity",values=qualitative_hcl(n = 6, h = c(0, 295), c = 80, l = 60, register = ) ) + labs(x=" ", y="CPUE (kg/hrs fishing) / number of traps pr area") + theme_bw() + 
  theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank())
t2.plot

tiff(file="./figs/diel_act.tiff", width=3000, height=3000, res=400, pointsize=10, compression=c("none"))
t2.plot
dev.off()
## quartz_off_screen 
##                 2

Gregariousness

t3.plot <- ggplot(subset(catch.traits, Gregariousness!="#N/A" & gear=="TB"), aes(Gregariousness, CPUEw/NArea))

t3.plot <- t3.plot + geom_col(aes(fill=Gregariousness)) + facet_grid(vars(id), vars(survey_m)) + scale_fill_manual(name="Gregariousness",values=qualitative_hcl(n = 6, h = c(0, 295), c = 80, l = 60, register = ) ) + labs(x=" ", y="CPUE (kg/hrs fishing) / number of traps pr area") + theme_bw() + 
  theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank())
t3.plot

tiff(file="./figs/gregariousness.tiff", width=3000, height=3000, res=400, pointsize=10, compression=c("none"))
t3.plot
dev.off()
## quartz_off_screen 
##                 2

Biodiversity (species densities, trophic diversity etc)

Aggregated Species density per area for (traps only)

Species per area pr hour of fishing, traps only.

The highest number of species was observed in area 2 (36), followed by area 3, 6 and 7 (25 in each).

Species density (number of species pr hour of fishing - soak time of traps) ranged from 0.008 (areas 1, 2 and 6) to 0.019 (area 4). Area 4 had the highest species density across all surveys (figure below), while the area with the lowest species density varied by survey: area 1 in November 2012 , area 7 in May 2013 , and area 2 in November 2013 .

sp.pr.area <- subset(catch.unmod, gear=="TB") %>% group_by(id, survey_m) %>% summarise(n_distinct(species), n_distinct(ss), sum(Fhrs))
colnames(sp.pr.area) <- c("id", "survey","Nsp","NTraps", "Fhrs")

sp.pr.area$Species_by_NoStations <- (sp.pr.area$Nsp-1) / sp.pr.area$NTraps
sp.pr.area$Species_Hrs <- (sp.pr.area$Nsp-1) / sp.pr.area$Fhrs
write.csv2(sp.pr.area, "./figs/species_density_area.csv")

#calculate species density pr area (all years combined)
sp.pr.area2 <- subset(catch.unmod, gear=="TB") %>% group_by(id) %>% summarise(n_distinct(species), n_distinct(ss), sum(Fhrs))
colnames(sp.pr.area2) <- c("id","Nsp","NTraps", "Fhrs")
sp.pr.area2$Species_by_NoStations <- (sp.pr.area2$Nsp-1) / sp.pr.area2$NTraps
sp.pr.area2$Species_Hrs <- (sp.pr.area2$Nsp-1) / sp.pr.area2$Fhrs
formattable(round(sp.pr.area2, digits=3))
id Nsp NTraps Fhrs Species_by_NoStations Species_Hrs
1 18 74 2036.710 0.230 0.008
2 36 192 4483.859 0.182 0.008
3 25 67 1471.675 0.358 0.016
4 12 34 593.399 0.324 0.019
5 19 94 2097.324 0.191 0.009
6 25 121 3114.607 0.198 0.008
7 25 92 2485.890 0.261 0.010 
#' Species Density plot
species.plot<-ggplot(sp.pr.area, aes(x=id, y=Species_Hrs)) + geom_bar(stat="identity", fill="steelblue") +theme_bw() + scale_x_reverse(breaks=c(1,2,3,4,5,6,7)) 
sphrs_unmod<- species.plot  + coord_flip() + theme(legend.title = element_text(size=14, face="bold")) + theme(legend.text = element_text(size=12)) + ylab("Number of species caught in traps / median hours of  fishing") + xlab("Management areas from N (1) to S (7)") + facet_wrap(~survey)
sphrs_unmod

ggsave("./figs/species_density.png")
## Saving 7 x 5 in image

Species density pr station - Zero Adjusted GAM model

The variability and relative difference in rank of species diversity per area as observed in the FIGURE ABOVE was confirmed by the zero-adjusted GAM model (AIC: -4161, df: 12) were all surveys were significant factors in the model. In addition, area 1, 2, 3 and 4 were also significant factors affecting species density, while the non-significant factors were: depth, and areas 5, 6 and 7.

spue.df <-  subset(catch.unmod, gear=="TB") %>% group_by(ss) %>% summarise(first(survey), first(id), n_distinct(species), first(depth), sum(Fhrs))
colnames(spue.df) <- c("id", "survey", "area", "Nsp", "depth", "Fhrs")
spue.df$Species_Hrs <- spue.df$Nsp / spue.df$Fhrs
#hist(spue.df$spue)

# zero-adjusted model:
sd.mod<-gamlss(Species_Hrs~depth+factor(area)+factor(survey),family=ZAGA, data=spue.df) # AIC: -4161
## GAMLSS-RS iteration 1: Global Deviance = -4185.753 
## GAMLSS-RS iteration 2: Global Deviance = -4185.753
# gamlss(Species_Hrs~factor(area)+factor(survey),family=ZAGA, data=spue.df)  # AIC: - 4163 - slightly better than the one above, but want to show the one with 'depth' as a factor. 
summary(sd.mod)
## ******************************************************************
## Family:  c("ZAGA", "Zero Adjusted GA") 
## 
## Call:  gamlss(formula = Species_Hrs ~ depth + factor(area) +  
##     factor(survey), family = ZAGA, data = spue.df) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                         Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)           -3.0456597  0.0282803 -107.695  < 2e-16 ***
## depth                  0.0003143  0.0003925    0.801 0.423536    
## factor(area)2          0.0952550  0.0257799    3.695 0.000238 ***
## factor(area)3          0.1083585  0.0319608    3.390 0.000740 ***
## factor(area)4          0.0834147  0.0390748    2.135 0.033148 *  
## factor(area)5          0.0540349  0.0296866    1.820 0.069184 .  
## factor(area)6         -0.0337542  0.0278053   -1.214 0.225200    
## factor(area)7         -0.0464743  0.0294473   -1.578 0.114993    
## factor(survey)2013002  0.2071281  0.0174699   11.856  < 2e-16 ***
## factor(survey)2013005  0.2470968  0.0190004   13.005  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.66955    0.02708  -61.66   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  logit 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)    -16.2      129.2  -0.125      0.9
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  674 
## Degrees of Freedom for the fit:  12
##       Residual Deg. of Freedom:  662 
##                       at cycle:  2 
##  
## Global Deviance:     -4185.753 
##             AIC:     -4161.753 
##             SBC:     -4107.594 
## ******************************************************************

Functional diversity

The ‘dbFD’ function only works when run in-line, not the whole code-chunk.

We used the same approach as Stuart-Smith et al to exclude very few observations of a single species. and excluded all species observations occuring in 2 or fewer stations (traps).

To avoid dbFD crashing ‘m’ was set in the ‘dbFD’ call (the number of PoCA axes kept as traits to do the FRic analysis). Several levels of ‘m’ were tested, and this could probably be tested further (to higher levels of ‘m’)

Results for n_spec > 2, m=10

Results based on the output from the dbFD analysis:

No area came out with the highest or lowest FD across all four metrics, but some patterns emerged:

  • Area 1 was highest in FEve and FDiv, but lowest in FDis and RaoQ
  • Area 2 was consistently high in all metrics, and highest in RaoQ
  • Area 3 was high inn all metrics, but never the highest
  • Area 4 was lowest in RaoQ, low in FDis, medium in FDiv and FEve
  • Area 5 was high in all metrics, but never the highest
  • Area 6 was high in FEve, but the lowest in all other three
  • Area 7 was lowest in FEve, low in FDis, medium in FDiv and RaoQ

RaoQ is the most used metric (the one used by Stuart-Smith et al), and this is the metric we will use and discuss in the MS. It supports the analysis of species densities for Area 2, although Area 6 has the lowest RaoQ which indicates that although the number of species is higher, their functional diversity is lower (or what we sampled was lower) compared to the other areas.

Stuart-Smith tests various RaoQs by taking out one trait at a time and this way ranking the traits according to R2.

#ss_list <- levels(catch1$survey)
sp_occur <- catch1  %>% group_by(Sci_name) %>% summarise(n())
colnames(sp_occur) <- c("Sci_name", "n_spec")

# wrangle traits into a traitstable that can be used in FD
missing <- setdiff( x=traits$x1, y=subset(sp_occur, n_spec>2)$Sci_name)
t5 <- traits[!traits$x1 %in% missing,]
tr <- t5[,c(4:10)]
rownames(tr) <- t5$x1
tr[(tr=="#N/A")] <- c(NA)
tr <- droplevels(tr)
tr$diel<-as.integer(tr$Diel.Activity)
tr$diel[(tr$diel==2)] <- c(0)
tr <- tr[,c(1,2,3,4,6,7,8)]

# wrangle 'catch1' to aggregate CPUE pr area (row) pr species (columns)
missing <- setdiff( x=catch1$Sci_name, y=rownames(tr))
c5 <- catch1[!catch1$Sci_name %in% missing,]
c_melt <- melt(c5[,c(28,11,22)], id=c("id", "Sci_name"))

#convert abundance to integer for better calculation performance
med_int1000 <- function(x) as.integer(median(x)*1000)
abund <- dcast(c_melt, id ~ Sci_name, value.var="value", fun=med_int1000) 
# to calculate with raw abundance data, use line below:
#abund <- dcast(c_melt, id ~ Sci_name, value.var="value", fun=median)
abund[is.na(abund)] <- 0
abund <- abund[,c(2:length(abund))]

#calculate FD using dbFD function, m set to 10 to avoid crach (int overload in QHULL)
FuncDiv <- dbFD(tr, abund, corr = c("lingoes"), m=10)
## Species x species distance matrix was not Euclidean. Lingoes correction was applied. 
## FRic: Dimensionality reduction was required. The last 49 PCoA axes (out of 59 in total) were removed. 
## FRic: Quality of the reduced-space representation (based on corrected distance matrix) = 0.4209541
FDt <- as.data.frame(rbind(FuncDiv$FEve, FuncDiv$FDiv, FuncDiv$FDis, FuncDiv$RaoQ))
rownames(FDt) <- c("FEve", "FDiv", "FDis", "RaoQ")
colnames(FDt) <- c("A1", "A2", "A3", "A4", "A5", "A6", "A7")

FDt <- as.data.frame(t(round(FDt, digits=3)))

# create a nice table of FD indicies for the 7 management areas
formattable(FDt, list(    FEve = color_tile("transparent", "lightpink"),     FDiv = color_tile("transparent", "lightpink"),     FDis = color_tile("transparent", "lightpink"),    RaoQ = color_tile("transparent", "lightpink")))
FEve FDiv FDis RaoQ
A1 0.682 0.905 0.189 0.051
A2 0.637 0.849 0.263 0.079
A3 0.644 0.710 0.239 0.074
A4 0.614 0.745 0.200 0.048
A5 0.627 0.785 0.247 0.073
A6 0.665 0.650 0.179 0.048
A7 0.516 0.713 0.197 0.064

Rerunning the FD model excluding one trait at a time

Need to reduce ‘m’ to 9 to avoid crash:

Zero distance(s)Error in convhulln(tr.FRic, “FA”) : Received error code 2 from qhull. Qhull error: QH6114 qhull precision error: initial simplex is not convex. Distance=1.4e-16

Removing the MaxLenght trait from the model increased R2 by 0.047 compared to the model with all traits included. R2 was reduced when removing any of the other six traits from the model.

Conclusion

  • Use the model without MaxLenght
  • Removing Maxlength makes sence because the data included a few catches of sharks with very large MaxLengt (>2m) that skews the MaxLenght distribution. Excluding MaxLength creates an analysis around traits that are more common in distribution across species.
  • Also, the sharks are rowing predators, not so site attached (but that also applies to some pelagic species, especially the large ones, so the same argument applies for these - excluding large body size reduces the importance of these non-site attached species)
# create data frame with original RaoQ and FRic
Raoq <- as.data.frame(c(FuncDiv$RaoQ, FuncDiv$qual.FRic))
colnames(Raoq) <- c("All")
rownames(Raoq) <- c("A1", "A2", "A3", "A4", "A5", "A6", "A7", "R2")


# select tr, but excluding one variable at a time
#tr[,c(1,!2,3,4,5,6,7)]
for (i in 1:7) { 
  FDout <- dbFD(tr[,c(setdiff(1:7, i))], abund, corr = c("lingoes"), m=9)
  raoq1 <- as.data.frame(c(FDout$RaoQ, FDout$qual.FRic))
  colnames(raoq1) <- colnames(tr)[i]
  Raoq <- cbind(Raoq, raoq1)
}
## Species x species distance matrix was not Euclidean. Lingoes correction was applied. 
## FRic: Dimensionality reduction was required. The last 37 PCoA axes (out of 46 in total) were removed. 
## FRic: Quality of the reduced-space representation (based on corrected distance matrix) = 0.4684034 
## Species x species distance matrix was not Euclidean. Lingoes correction was applied. 
## FRic: Dimensionality reduction was required. The last 42 PCoA axes (out of 51 in total) were removed. 
## FRic: Quality of the reduced-space representation (based on corrected distance matrix) = 0.3984771 
## Species x species distance matrix was not Euclidean. Lingoes correction was applied. 
## FRic: Dimensionality reduction was required. The last 50 PCoA axes (out of 59 in total) were removed. 
## FRic: Quality of the reduced-space representation (based on corrected distance matrix) = 0.3663967 
## Species x species distance matrix was not Euclidean. Lingoes correction was applied. 
## FRic: Dimensionality reduction was required. The last 50 PCoA axes (out of 59 in total) were removed. 
## FRic: Quality of the reduced-space representation (based on corrected distance matrix) = 0.3612455 
## Species x species distance matrix was not Euclidean. Lingoes correction was applied. 
## FRic: Dimensionality reduction was required. The last 48 PCoA axes (out of 57 in total) were removed. 
## FRic: Quality of the reduced-space representation (based on corrected distance matrix) = 0.4029848 
## Species x species distance matrix was not Euclidean. Lingoes correction was applied. 
## FRic: Dimensionality reduction was required. The last 50 PCoA axes (out of 59 in total) were removed. 
## FRic: Quality of the reduced-space representation (based on corrected distance matrix) = 0.3538643 
## Species x species distance matrix was not Euclidean. Lingoes correction was applied. 
## FRic: Dimensionality reduction was required. The last 50 PCoA axes (out of 59 in total) were removed. 
## FRic: Quality of the reduced-space representation (based on corrected distance matrix) = 0.4132985
formattable(round(Raoq[8,], digits=3))
All MaxLength Trophic.Level Trophic.group Water.column Habitat Gregariousness diel
R2 0.421 0.468 0.398 0.366 0.361 0.403 0.354 0.413 
formattable(round(Raoq[c(1:7),], digits=3), list(    MaxLength = color_tile("transparent", "lightpink"),     Trophic.Level = color_tile("transparent", "lightpink"),     Trophic.group = color_tile("transparent", "lightpink"),    Water.column = color_tile("transparent", "lightpink"), Habitat = color_tile("transparent", "lightpink"),    Gregariousness = color_tile("transparent", "lightpink"), diel = color_tile("transparent", "lightpink")))
All MaxLength Trophic.Level Trophic.group Water.column Habitat Gregariousness diel
A1 0.051 0.055 0.065 0.049 0.040 0.063 0.040 0.059
A2 0.079 0.102 0.103 0.059 0.059 0.093 0.070 0.089
A3 0.074 0.078 0.092 0.061 0.061 0.100 0.062 0.077
A4 0.048 0.063 0.063 0.050 0.031 0.046 0.039 0.057
A5 0.073 0.085 0.093 0.060 0.055 0.091 0.061 0.082
A6 0.048 0.050 0.058 0.036 0.042 0.064 0.046 0.051
A7 0.064 0.076 0.083 0.057 0.050 0.067 0.060 0.070 
write.csv2(Raoq, "./figs/Raoq_table.csv")

Results table

This table combines CPUE,

st2 <- catch %>% group_by(ss) %>% summarise(first(id), first(Fhrs), first(gear), sum(weight))
colnames(st2) <- c("ss", "id", "Fhrs","gear", "weight")
st2$CPUEw <- st2$weight/st2$Fhrs

results.table <- data.frame(area=(1:7))
results.table <- cbind(results.table,(subset(st2, gear=="TB") %>% group_by(id) %>% summarise( mean(CPUEw)))[,2])
results.table <- cbind(results.table,(subset(st2, gear=="GN") %>% group_by(id) %>% summarise( mean(CPUEw)))[,2])
results.table <- cbind(results.table,(subset(catch.traits, gear=="TB") %>% group_by(id) %>% summarise( mean(TrophicLevel, na.rm = TRUE)))[,2])
results.table <- cbind(results.table,(subset(catch.traits, gear=="GN") %>% group_by(id) %>% summarise( mean(TrophicLevel, na.rm = TRUE)))[,2])
results.table<- cbind(results.table,(subset(catch.traits, gear=="TB") %>% group_by(id) %>% summarise( mean(as.integer(Gregariousness), na.rm = TRUE)))[,2])
results.table<- cbind(results.table,(subset(catch.traits, gear=="GN") %>% group_by(id) %>% summarise( mean(as.integer(Gregariousness), na.rm = TRUE)))[,2])
colnames(results.table) <- c("area", "TrapCPUE", "GillnetCPUE", "TrophicL_Traps", "TrophicL_Gillnets", "Greg_Traps", "Greg_Gillnets" )

TMean <- subset(catch.traits, gear=="TB" & TrophicLevel!="NA") %>% group_by(id, survey_m) %>% summarise(TM=mean(na.omit(TrophicLevel)))

results.table$SpeciesDensity <- sp.pr.area2$Species_Hrs
results.table$RaoQ <- Raoq$MaxLength[1:7]

formattable(round(results.table, digits = 3))
area TrapCPUE GillnetCPUE TrophicL_Traps TrophicL_Gillnets Greg_Traps Greg_Gillnets SpeciesDensity RaoQ
1 0.117 1.639 3.825 3.933 1.510 2.412 0.008 0.055
2 0.149 1.176 3.815 3.891 1.478 2.261 0.008 0.102
3 0.136 0.309 3.844 4.140 1.553 2.500 0.016 0.078
4 0.059 0.698 3.765 4.284 1.357 2.107 0.019 0.063
5 0.079 0.726 3.797 3.853 1.627 2.176 0.009 0.085
6 0.135 0.270 3.762 3.927 1.531 2.286 0.008 0.050
7 0.098 0.683 3.828 3.612 1.565 1.900 0.010 0.076 
write.csv2(round(results.table, digits = 3), "./figs/results_table.csv")

Statistical analysis CPUE catches

Testing for normality of CPUE catches

Shapiro test for normality and Q-Q plots shows that data is non-normal and zero-inflated.

catch.df<-subset(catch1, number>0 & weight>0)

#' test for normality
shapiro.test(catch.df$CPUEw)
## 
##  Shapiro-Wilk normality test
## 
## data:  catch.df$CPUEw
## W = 0.5481, p-value < 2.2e-16
shapiro.test(catch.df$CPUEn)
## 
##  Shapiro-Wilk normality test
## 
## data:  catch.df$CPUEn
## W = 0.38853, p-value < 2.2e-16
#' Q-Q plots
library("qqplotr")
## 
## Attaching package: 'qqplotr'
## The following objects are masked from 'package:ggplot2':
## 
##     stat_qq_line, StatQqLine
gg <- ggplot(data = catch.df, mapping = aes(sample = CPUEw,  fill = factor(id))) + stat_qq_band() + stat_qq_line() + stat_qq_point(size=0.2) + facet_wrap(~ id, scales="free") +  labs(x = "CPUEw predicted", y = "CPUEw Observed") + ggtitle("Q-Q plots CPUEw")
gg

### Test for difference in CPUE between areas (non-parametric test) Both CPUEw and CPUEn are significantly different between the areas, but only CPUEn in the pairwise tests (nemenyi post-hoc) for areas 4:3 and 6:4. THis points to that Area 4 is the ‘odd one out’ in terms of CPUE along the Sudanses coast.

Significant difference for CPUEW between May 2013 survey and the Nov 2012 and Nov 2013 surveys, but not between the Nov12 and Nov13 surveys

#### NOT RUN BECAUSE ZERO-ADJUSTED GAM MODEL IS BETTER!

#' 1 - test for difference between areas
kruskal.test(catch.df$CPUEw, catch.df$id)
kruskal.test(catch.df$CPUEn, catch.df$id)

#' Post-Hoc tests of Kruskal - Wallis results
library("PMCMR", lib.loc="/Library/Frameworks/R.framework/Versions/3.5/Resources/library")
posthoc.kruskal.nemenyi.test(catch.df$CPUEw, catch.df$id, "Chisq")
posthoc.kruskal.nemenyi.test(catch.df$CPUEn, catch.df$id, "Chisq")
#' significant difference in CPUEn between area 6 and area 4

#' 2 test for differences between surveys
kruskal.test(catch.df$CPUEw, catch.df$survey)
kruskal.test(catch.df$CPUEn, catch.df$survey)

posthoc.kruskal.nemenyi.test(catch.df$CPUEw, catch.df$survey, "Chisq")
posthoc.kruskal.nemenyi.test(catch.df$CPUEn, catch.df$survey, "Chisq")

Zero-inflated model (better approach)

This method is more appropriate than the Kruskal-Wallis approach above as it evaluates the effects of all factors at the same time, instead of evaluating the effects of area separate from survey.

First, developed a model for CPUEw including depth, area, survey and gear, which performed better than a simpler model excluding gear.

Based on reviewer comments included survey as a random factor (using two different approaches: ´random´and ´re´), although this excludes evaluation of the potential significant effects of the surveys on the results.

## GAMLSS-RS iteration 1: Global Deviance = 223.3244 
## GAMLSS-RS iteration 2: Global Deviance = 223.3244
## GAMLSS-RS iteration 1: Global Deviance = 226.7143 
## GAMLSS-RS iteration 2: Global Deviance = 226.8294 
## GAMLSS-RS iteration 3: Global Deviance = 226.83
## GAMLSS-RS iteration 1: Global Deviance = 226.7306 
## GAMLSS-RS iteration 2: Global Deviance = 226.8359 
## GAMLSS-RS iteration 3: Global Deviance = 226.8359
##            df      AIC
## mod2 11.00294 248.8359
## mod3 11.00002 248.8359
## mod  13.00000 249.3244
## ******************************************************************
## Family:  c("ZAGA", "Zero Adjusted GA") 
## 
## Call:  gamlss(formula = CPUEw ~ depth + factor(id) + factor(gear) +  
##     random(factor(survey)), family = ZAGA, data = c2) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -1.409137   0.152478  -9.242  < 2e-16 ***
## depth           0.004225   0.002393   1.765 0.077824 .  
## factor(id)2    -0.026764   0.149702  -0.179 0.858143    
## factor(id)3    -0.008394   0.201729  -0.042 0.966818    
## factor(id)4    -0.045095   0.206644  -0.218 0.827291    
## factor(id)5    -0.444315   0.168467  -2.637 0.008468 ** 
## factor(id)6    -0.343869   0.166375  -2.067 0.038976 *  
## factor(id)7    -0.531138   0.164365  -3.231 0.001267 ** 
## factor(gear)TB -0.385532   0.100217  -3.847 0.000126 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.09767    0.02207   4.425 1.06e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  logit 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.70099    0.06266  -11.19   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  1149 
## Degrees of Freedom for the fit:  11.00294
##       Residual Deg. of Freedom:  1137.997 
##                       at cycle:  3 
##  
## Global Deviance:     226.83 
##             AIC:     248.8359 
##             SBC:     304.3639 
## ******************************************************************

## ******************************************************************
##   Summary of the Randomised Quantile Residuals
##                            mean   =  -0.03174958 
##                        variance   =  1.047609 
##                coef. of skewness  =  0.3344842 
##                coef. of kurtosis  =  3.863233 
## Filliben correlation coefficient  =  0.9920283 
## ******************************************************************

## GAMLSS-RS iteration 1: Global Deviance = 42.87 
## GAMLSS-RS iteration 2: Global Deviance = 42.9993 
## GAMLSS-RS iteration 3: Global Deviance = 43.0024 
## GAMLSS-RS iteration 4: Global Deviance = 43.0025
## ******************************************************************
## Family:  c("ZAGA", "Zero Adjusted GA") 
## 
## Call:  gamlss(formula = CPUEw ~ depth + factor(id) + factor(gear) +  
##     random(factor(survey)), nu.formula = ~depth + factor(id) +  
##     factor(gear) + random(factor(survey)), family = ZAGA,      data = c2) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -1.409401   0.152483  -9.243  < 2e-16 ***
## depth           0.004232   0.002394   1.768 0.077343 .  
## factor(id)2    -0.026754   0.149706  -0.179 0.858199    
## factor(id)3    -0.008466   0.201730  -0.042 0.966531    
## factor(id)4    -0.044849   0.206657  -0.217 0.828231    
## factor(id)5    -0.444299   0.168469  -2.637 0.008472 ** 
## factor(id)6    -0.343919   0.166378  -2.067 0.038953 *  
## factor(id)7    -0.531118   0.164369  -3.231 0.001268 ** 
## factor(gear)TB -0.385442   0.100216  -3.846 0.000127 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.09767    0.02207   4.425 1.06e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  logit 
## Nu Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -3.157336   0.341870  -9.235   <2e-16 ***
## depth          -0.002249   0.003735  -0.602    0.547    
## factor(id)2     0.058979   0.244381   0.241    0.809    
## factor(id)3     0.268465   0.305364   0.879    0.379    
## factor(id)4     0.546133   0.364997   1.496    0.135    
## factor(id)5     0.426533   0.280090   1.523    0.128    
## factor(id)6    -0.296271   0.268989  -1.101    0.271    
## factor(id)7     0.110213   0.273484   0.403    0.687    
## factor(gear)TB  2.833560   0.286698   9.883   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  1149 
## Degrees of Freedom for the fit:  21.75223
##       Residual Deg. of Freedom:  1127.248 
##                       at cycle:  4 
##  
## Global Deviance:     43.00247 
##             AIC:     86.50693 
##             SBC:     196.2828 
## ******************************************************************
## GAMLSS-RS iteration 1: Global Deviance = 58063.27 
## GAMLSS-RS iteration 2: Global Deviance = 58063.27
## Warning in summary.gamlss(mod5): summary: vcov has failed, option qr is used instead
## ******************************************************************
## Family:  c("ZAIG", "Zero adjusted IG") 
## 
## Call:  gamlss(formula = CPUEw ~ depth + factor(id) + factor(gear) +  
##     random(factor(survey)), family = ZAIG, data = c2) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)
## (Intercept)    -1.451e+00  1.158e+04       0        1
## depth           3.467e-03  1.612e+02       0        1
## factor(id)2    -3.445e-02  1.097e+04       0        1
## factor(id)3    -3.955e-02  1.425e+04       0        1
## factor(id)4    -4.657e-02  1.519e+04       0        1
## factor(id)5    -4.383e-01  1.232e+04       0        1
## factor(id)6    -3.308e-01  1.200e+04       0        1
## factor(id)7    -5.247e-01  1.209e+04       0        1
## factor(gear)TB -3.009e-01  8.080e+03       0        1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##              Estimate Std. Error   t value Pr(>|t|)    
## (Intercept) 3.969e+01  2.950e-07 134552198   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  logit 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.70099    0.06266  -11.19   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms may not be reliable. 
## ------------------------------------------------------------------
## No. of observations in the fit:  1149 
## Degrees of Freedom for the fit:  11.82643
##       Residual Deg. of Freedom:  1137.174 
##                       at cycle:  2 
##  
## Global Deviance:     58063.27 
##             AIC:     58086.92 
##             SBC:     58146.61 
## ******************************************************************

Significant factors for both models with random effects (mod2 and mod3): * Depth * Area (1, 5, 6 and 7) * Gear (traps)